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ABSTRACT 

We present a series of numerical sunspot models addressing the subsurface field and flow structure in 
up to 16 Mm deep domains covering up to 2 days of temporal evolution. Changes in the photospheric 
appearance of the sunspots are driven by subsurface flows in several Mm depth. Most of magnetic 
field is pushed into a downflow vertex of the subsurface convection pattern, while some fraction of the 
flux separates from the main trunk of the spot. Flux separation in deeper layers is accompanied in the 
photosphere with light bridge formation in the early stages and formation of pores separating from 
the spot at later stages. Over a time scale of less than a day we see the development of a large scale 
flow pattern surrounding the sunspots, which is dominated by a radial outflow reaching about 50% of 
the convective rms velocity in amplitude. Several components of the large scale flow are found to be 
independent from the presence of a penumbra and the associated Evershed flow. While the simulated 
sunspots lead to blockage of heat flux in the near surface layers, we do not see compelling evidence 
for a brightness enhancement in their periphery. We further demonstrate that the influence of the 
bottom boundary condition on the stability and long-term evolution of the sunspot is significantly 
reduced in a 16 Mm deep domain compared to the shallower domains considered previously. 
Subject headings: Sun: surface magnetism - sunspots - MHD - convection 



1. INTRODUCTION 

The subsurface structure of sunspots has been of sub- 
ject of theoretical investigations for several decades. The 
two possible extrem e s of m agnetic configurations were 
discussed by Parker (19791: a monolithic configuration 
vs. a clusters of individual flux tubes that is kept to- 
gether by converging flows in a suitable depth. On the 
observational side evidence is inconclusive. Direct he- 
lioseismic measurements of the magnetic field structure 
with an accuracy to determine the differences between 
monolithic and cluster models are currently out of reach. 
More promising are measurements of the subsurface flow 
structure, but also th ere results are inconclus ive. Time 
distance inversions by |Zhao et al.| ( |2001[ |2010|) point to- 
ward inflows around sunspots m an intermediate depth 
range from about 1.5 to 5 Mm (and corresponding down- 
flows underneath sunspots), which would be consistent 
with the expectations from a cluster model. On the other 



hand r ecent inversions presented by Gizon et al, 
2010b ) show outflows in the upper most 4.5 Mm 



(2009 



At photospheric levels most sunspots are surrounded 
by large scale outflows (called "moat flows" ) with ampli- 
tudes of a few 100 ms _1 that were first found through 



tracking of m agnetic features (Sheeley |1969[ |Harvey fc 



Harvey 1973), Doppler measurements (Shceley 19721 and 



later fielioseismic measurements ( jGizonet al.|2000[ ). Sev- 
eral recent investigations focused on possible conne ctions 
between the Evershed flow and moat flow (see e.g. Sainz 
Dalda & Martinez Pillet|2005l|Cabrera Solana et al. 2006 



Var gas Dommguez et al.||20081 |Zuccarello et al 



2009 



Vargas Dominguez et al. 2010p , so far the observationa 



evidence is not clear enough to either proof or disproof a 
connection. 

The subsurface field and flow structure has been also 
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addressed through mea ns of numerical models. 2D ax 
isymmetric models by Hurlburt fc Rucklidge| |2000), 
and 



Botha et al. (2006) 



Botha et al. 



scale flow patterns around sunspots 



2008P showlarge 
.'he typical result 



is a converging collar flow in the proximity of the spot 
and a diverging flow further out, a situations which is 
essentially in agree ment with the cluster mo del as well 



as the inversions b y Zhao et al. (2001 2010). Recently 



Botha et al. (2011 ) expanded this work to 3D and found 



comparable results with regard to the axisymmetric flow 
components. Differences occurred with regard to the pro- 
cess of sunspot decay: the azimuthal cell structure allows 
for flux to escape from the central flux concentration even 
if the average flow is converging. At this point neither 
the 2D axisymmetric nor the 3D simulations contain a 
penumbra and the connection of the larger scale deep 
seated flows to photospheric flows remains an open ques- 
tion. 

Over the past five years there has been a substantial 
progress in 3D MHD models that include a realistic equa- 
tion of state and radiative transfer. Owing to the wide 
range of length and time scales involved in sunspot struc- 
ture and evolution it is currently not possible to address 
all aspects of sunspot structure and evolution in a single 
numerical simulation. The formation, evolution and de- 
cay of pore-si z e flux c oncentrations has been modeled by 



Bercik et al. (2003); Cameron et al. (2007); Kitiashvili 



et al. "(2010), typically resulting in converging flows in 
the proximity of the pore. The focus of re c ent numeri 
cal models such as|Schussler fc Vogler (2006), Heinemann 



|et^([2Tj0T| , |RempeTet al.| ( |2009b[ ),and |Kitiashvili et al 
( |2009| was primarily the sunspot fine structure and ori- 
gin of the Evershed flow. To this end those models fo- 
cused on smaller subsections and rather short t emporal 
evolution of a few hours. The sim ulations by |Rempel| 



et al. (2009a); Rempel (2010 2011) were the first MHD 
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simulations with a sufficient domain size in the horizontal 
direction to capture complete sunspots, while still resolv- 
ing sunspot fine structure. However, the vertical domain 
size of 6.f 44 Mm remained still too shallow and the over- 
all time evolution of up to 6 hours too short to properly 
address the subsurface structure of the sunspot as well 
as development of moat flows. While these simulations 
present a unified picture for the magneto convective ori- 
gin of sunspot fine structure they leave several funda- 
mental aspects with regard to the subsurface structure 
unanswered. 

Evidence for a larger scale diverging flow is ahead ; 
present in the simula t ions of i Heinemann et al. ( 2007[ 
Rempel et al] ( |2009b[ ) , |Rempel et al.| ( |2UU9a[ ), and |Rem 
pell2Ull ), but the overall development of this flow pat 
tern was heavily influenced by either the domain size or 
the rather short time span covered. We will relax most 
of these constraints in this investigation by focusing on 
numerical simulations with lower resolution, but deeper 
domains (up to 16 Mm) and much longer temporal evo- 
lution (up to 2 days). Simulations in deeper domains are 
less influenced by the bottom boundary condition and 
allow us to study processes related to the fragmentation 
of sunspots and development of large scale flow patterns. 

In Sect. [2] we describe the setup of the numerical sim- 
ulations used in this investigation. Sect. [3] presents the 
time evolution and structure of the subsurface magnetic 
field and flow structure. The influence of the (unfortu- 
nately unavoidable) bottom boundary condition is ana- 
lyzed in Sect. |4j the findings are discussed in Sect. [5]and 
summarized in Sect. |6] 

2. NUMERICAL MODELS 

The details of the n umerical model an d un derlying 
physics are d escribed in |Vogler et al. (2005) and Rempel 



et all fl2009 bl . Models pre s ented h ere use a setup similar 



to those of Rempel et al.| ( |2009a[ ), but differ in domain 
size, resolution as well as boundary conditions. 

We report here on a series of numerical simulations 
which were performed in a computational domain with 
a horizontal extent of 49.152 Mm and a vertical extent 
of 16.384 Mm. We used a rather low resolution of 96 km 
horizontally and 32 km vertically to allow for long sim- 
ulation runs covering up to 2 days of solar time. The 
simulations were started from a non-magnetic convec- 
tion simulation, which was evolved for 42 hours to al- 
low for thermal relaxation of the stratification. We then 
inserted into the domain an axisymmetric, self-similar 
sunspot with 1.2 ■ 10 22 Mx flux and a field strength of 
16 kG at the bottom of the domain. While the domain is 
periodic in the horizontal direction, the top boundary is 
closed for the vertical mass flux and slip free for horizon- 
tal motions. The magnetic field is matched to a potential 
field extrapolation. At the bottom boundary the mag- 
netic field is vertical and we consider here two formula- 
tions with regard to flows. The first boundary condition 
sets all velocity components to zero in regions where the 
field strength exceeds 5 kG, while it allows for convective 
motions to cross the bottom boundary everywhere else: 
in outflow regions all velocity gradients are set to zero 
and thermodynamic variables are extrapolated into the 
ghost cells while inflows have purely vertical velocities 
and a prescribed entropy to maintain the solar energy 
flux. The second boundary condition is an open bound- 



ary everywhere in the domain regardless of the magnetic 
field strength. In contrast to the first boundary condi- 
tion we set the gradient of all mass flux components to 
zero, allowing also for horizontal flows in inflow regions 
(a zero gradient of the horizontal flow velocity instead of 
mass flux turns out to be unstable on long time scales). 
While the first boundary condition strongly suppresses 
the decay of the sunspot by limiting the motions at the 
foot point, the second boundary allows for horizontal ex- 
change between magnetized and unmagnetized regions 
as well as vertical flows changing the mass content of the 
sunspot. Note that the first boundary condition does 
not allow for heat flux entering the domain in strong 
field regions, while in the second case heat exchange is 
not prohibited by the boundary, although in general still 
strongly suppressed by the magnetic field. In both cases 
the energy radiated away in the sunspot umbra comes 
mostly from the heat content of the stratification, i.e., the 
umbra and stratification underneath slowly cool down as 
the simulation progresses. The most dramatic temper- 
ature changes happen during the first few hours of the 
simulations, later most of the energy comes from several 
Mm deep layers, where the heat capacity is sufficiently 
large to prevent strong temperature changes. A similar 
behavior is also fo und in simpler models such as Schiisslcr 



& Rempel {2005) 



in the following discussion we focus first on the simu- 
lation with the closed boundary condition in strong field 
regions, which leads to an overall more stable sunspot. 
We contrast these results with simulations using the sec- 
ond boundary condition in Sect. [4] 

We present in addition results from a simulation in a 
73.728 Mm wide and 9.216 Mm deep domain at a reso- 
lution of 48 km horizontally and 24 km vertically, which 
was evolved for a total of 24 hours. In this simulation we 
use a different top boundary condition, which enforces a 
more inclined magnetic field than the previous cases us- 
ing a potential field. To this end we enhanced the hori- 
zontal field components by a factor of 2 compared to a po- 
tential field extrapolation. Together with the higher res- 
olution this boundary leads to a stable extended penum- 
bra with average Evershed flows of about 4kms~ . The 
bottom boundary is closed in the inner most 8 Mm to pre- 
vent decay of the sunspot. The total flux of this sunspot 
is again 1.2 • 10 22 Mx, the field strength at the bottom 
boundary is 10 kG. We will present results from this sim- 
ulation with regard to axisymmetric mean flows and com- 
pare them to the previously mentioned simulation runs to 
investigate the influence of penumbra and Evershed flow 
on large scale flows surrounding sunspots. For a more 
detailed description of this simulation run and a series 
of simulations with the different top boundary condition 
we refer to a separate forthcoming publication. 

3. RESULTS 

3.1. Temporal evolution of subsurface magnetic field 

Figure [l] highlights the connection between the tempo- 
ral evolution of the subsurface magnetic field and appear- 
ance of the sunspot in the photosphere. The correspond- 
ing subsurface flow structure is presented in Figure [2] and 
discussed further in section [3~2l The left and middle col- 
umn show the magnetic field strength at a depth of 11.5 
and 3.3 Mm beneath the r = 1 level of the plage region 
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Figure 1. Temporal evolution of subsurface magnetic field strength at two depth levels (left column: z = —11.5 Mm, middle column: 
z = —3.3 Mm) and surface appearance of sunspot (right column: intensity). Top to bottom snapshots at t = 8.8 hours, t = 27.3 hours 
and t = 46.4 hours are shown. The magnetic field is swept into a downflow vertex of the subsurface convection pattern. A fraction of the 
initial magnetic flux is separating from the main trunk along downflow lanes. At the surface the separation of flux is accompanied with 
the formation of light bridges and pores surrounding the spot at later stages. An animation is available in the online material. 



surrounding the sunspot. The right column presents the 
surface intensity (gray intensity for vertical direction). 
We show top to bottom 3 time steps about 18 - 19 hours 
apart. The sunspots starts as a circular sunspot due 
to our initial condition, but is changing its shape sub- 
stantially throughout the simulated time span of about 
2 days. Several Mm beneath the photosphere most of 
the magnetic flux is collected into a downflow vertex (in 
about 10 Mm depth the intrinsic scale of convection be- 
comes larger than the diameter of the sunspot at that 
depth), while some of the magnetic flux separates from 



the main trunk along downflow lanes connecting to the 
vertex. A qualitatively similar picture of sunsp ot de- 
cay was also found in the recent 3D simulation of Botha 
et al. (2011). The flux separation driven by convective 



motions m deeper layers becomes manifest in the pho- 
tosphere as light bridges in the earlier stages and pores 
separating from the main spot during the later stages. 
While light bridges have a photospheric appearance very 
similar to penumbral filaments, their orig in is of fun- 
damentally d if ferent nature. As sho wn b y ISchussler fc 
|Vogler| ( |2006[ ) , |Rempel et al. ( 2009b ) and Rempel et al. 
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Figure 2. Evolution of large scale flows shown for the same snapshots presented in Figure [T] The left column present vertical velocity 
at a depth of z = —3.3 Mm, blue colors indicate upflows. The middle column presents radialflow velocity at a depth of z = —3.3 Mm, 
red colors indicate outflows. The right column presents azimuthal averages of the radial flow velocity, normalized by the convective rms 
velocity. The dotted horizontal line indicates the depth level of z = —3.3 Mm used in the left and middle column. Convection cells arrange 
in a ring-like fashion around the sunspot. As a consequence horizontal convective flows lead to a mean flow reaching about 50% of the 
convective rms velocity during the later stages of the simulation. 



(2009a I umbral dots and penumbral filaments originate 



trom magneto-convection in strong magnetic field. The 
observed reduction of field strength in these features is a 
consequence of overturning convection and does not re- 
quire the intrusion of "field free" plasma from beneath 
or outside the sunspot. The magneto-convective motions 
responsible for these structures are concentrated in the 
upper most 500 km beneath the photosphere. In con- 
trast to this light-bridges are the consequence of almost 
field free plasma entering from beneath and the associ- 
ated structure is deep reaching (several Mm). A good 
example is the snapshot at t = 27.3 hours shown in Fig- 
ure [T] The light-bridge entering the sunspot from the left 



side is clearly visible as an intrusion of field free plasma 
in 3.3 Mm depth and even the field in 11.5 Mm depth 
shows a similar signature. While narrow light bridges 
show the formation of a dark lane in their center (above 
the central upflow), wider light bridges break down into 
several granulation like convection cells. Good examples 
for the latter are present starting from t = 39 hours in 
the animation provided for Figure [T] in the online ma- 
terial. Light bridges with similar properties w ere also 



found in the simulations of Cheung et al. (2010). 



As a byproduct of the sunspot decay we find several 
pores surrounding the sunspot. The pores show a wide 
spread in life times and we see a very strong indication 
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Figure 3. Radial and vertical flow velocities at 3 depth levels for the snapshot at t = 46.4 hours. The top row shows radial, the bottom 
row vertical flow velocity. Left to right the depth levels z = —11.5 Mm, z = —3.3 Mm, and z = —1.3 Mm. The ring-like arrangement 
of convection cells is present at all depth levels, the diameter of the region in which the presence of the sunspot modifies the convection 
pattern is increasing with depth as the intrinsic scale of convection is increasing with pressure scale height. An animation is available in 
the online material. 



that the depth extent of the subsurface magnetic field 
structure is the most relevant factor determining their 
temporal evolution. Pores such as those found in the 
lower right corner of the snapshot at t = 46.4 hours in 
Figure [Rare connected to layers deeper than 10 Mm and 
are visible in the photosphere for almost 24 hours. 

3.2. Temporal evolution of subsurface velocity field 



penumbra and near surface flows are converging toward 
the spot with a velocity of several kins -1 . A similar be- 
havior is observed near pores, where the granules at the 
periphery of the umbra lead to a converging flow pattern 



(Wang fe Zirin 1992 Sobotka et al. 19991. A converg- 



in g how was also f ound in JVLHU simulations of a pore 

by " 



Cameron et al. 



Recently Rempel et al. (2009a) presented a numeri- 
cal simulation including an extended penumbra and Ev- 
ershed flows with average flow amplitudes faster than 
4kms _1 . Requirements for obtaining these results where 
a sufficiently high resolution to resolve penumbral fila- 
ments and the proximity of an opposite polarity spot 
leading to more horizontal magnetic field. In the simula- 
tions we discuss here we focus on an individual sunspot 
and our horizontal grid resolution is 96 km instead of 
32 km. Under these circumstances we do not obtain a 



( 2007 ) . We focus in the following dis 
cussion on large scale hows in the deeper domain, which 
develop on time scales of several hours to days. 

Figure [2] presents the time evolution of the subsurface 
flow fieldfor the same time steps shown in Figure [Tl 
The two columns on the left show vertical and radial 
velocity (with respect to the approximate center of the 
spot) in a depth of —3.3 Mm, which corresponds to the 
magnetic held evolution shown in the middle column of 
Figure [T] Here, and throughout this paper, upflows are 
represented by positive values and blue colors, outflows 
are represented by positive values and red colors. The 
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Figure 4. Properties of the azimuthally and temporally (15 hours from t = 35 hours to t = 50 hours) averaged flow in the simulation 
run without penumbra. Displayed are a) radial and b) vertical flow velocity relative to the convective rms velocity (positive values indicate 
out/up-flows). Panel c) presents the fraction of the mass flux present in the azimuthal mean component, panel d) the temperature 
perturbation relative to the mean stratification computed outside R = 25 Mm (i.e. corners of the domain) in units of the rms fluctuation 
(outside R = 25 Mm). It is evident that the main upflow coincides with an annulus of enhanced temperature between R = 8 and R = 16 Mm. 
The thermal signature of the large flow is comparable to that of convective flows. We indicated the field lines of the azimuthally averaged 
magnetic field for reference, the almost horizontal line near the top indicates the r = 0.1 level. 



right column presents azimuthal averages of the radial 
flow velocity as function of distance from the spot center 
(horizontal axis) and depth (vertical axis). Since convec- 
tive flow velocities vary quite substantially in a 16 Mm 
deep domain from a few 100 ms -1 near the base to a 
few kms -1 in the photosphere, we normalize here the 
azimuthal flow velocity by the convective rms velocity at 
the corresponding depth (averaged over the region out- 
side the sunspot). Over time we see the development of a 
radial mean flow with a weak inflow close to the sunspot 
and a large scale outflow in the region R > 10 Mm. To- 
ward the end of the simulation run the azimuthally aver- 
aged outflow reaches an amplitude of about 50% of the 
convective rms velocity, while the inflow reaches values 
around 15% (faster flows are found in the photosphere). 
Furthermore we see a strong trend of declining inflow and 
increasing outflow amplitude over time. The presence of 
the sunspot leads over time to a ring-like arrangement of 
convection cells in the periphery of the spot, w hich was 
already found in the m o re sha llow simulation of [Rempel] 
et al.| ( |2009a|); | Remp elj (|2 01 11) and the idealized simula- 
tions or |ljotha et al.| ( |20TT ). The result is a ring of upflow- 
ing material with overall reduced vertical flow velocity 



and fewer downflow lanes relative to convection further 
away from the sunspot (left column). The radial flow 
patterns (middle and right column) shows a converging 
flow in the proximity of the spot and a diverging flow fur- 
ther out. At earlier time steps the large scale flows near 
the bottom of the domain are less pronounced since the 
re-arrangement of the convection pattern requires several 
turn-over time scales and the latter increases substan- 
tially with depth. For example, the quantity H p /v TUls is 
about 1 hour in the middle of the domain, but 6-8 hours 
near the bottom of the domain. 

Figure [3] shows for a single time step (t = 46.4 hours) 
vertical and radial flow velocity for the depth levels of 
— 11.5, —3.3 and —1.3 Mm. The ring-like arrangement 
of convection cells is present at all depth levels shown, 
the scale of the resulting large scale flow patterns is in- 
creasing with depth as the intrinsic scale of convection is 
increasing proportional to the pressure scale height. This 
tendency is also manifest in the azimuthal average of the 
radial flow velocity toward the end of the simulation run 
(see Figure [2] for t = 46.4 hours). 

To analyze the large scale flows in more detail we fo- 
cus now on long-term time averages of the azimuthally 
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Figure 5. Same quantities as in Figure [4] for the simulation with penumbra and Evershed flow (see Figure [To] for an intensity snapshot 
and magnetogram from this simulation run). Differences compared to the previous case are restricted mostly to photosphcric layers (see 
indicated r = 0.1 level). While the previous case shows an inflow for R < 10 Mm and an outflow further out, radial flows are outward 
directed for all distances in the presence of a penumbra. The large scale outward directed flow with amplitudes around 0.5f rm s and the 
associated thermal perturbation is very similar to the previous case. Overall this points toward a large degree of independence of this flow 
component from the presence of a penumbra and Evershed flow. 



averaged radial and vertical flow velocity, mass flux and 
temperature perturbations. To this end we averaged 19 
snapshots between t — 35 and t = 50 hours (about 
50 minute cadence). The result is presented in Figure 
HI Velocities are normalized by the convective rms ve- 
locity, the temperature fluctuation is defined relative to 
the mean stratification, both, rms velocity and reference 
stratification are computed from the region surrounding 
the sunspot near the edges of the computational do- 
main. We indicate the extent of the sunspot through 
field lines of the azimuthally averaged field. The outer- 
most field line encloses 6.5 • 10 21 Mx flux. We contrast 
this result with a sunspot that has an extended penum- 
bra and Evershed flows reaching more than 4kms~ 1 on 
average in Figure [5] (see description in Sect. 2 for de- 
tails). All quantities are the same as in Figure 4" except 
for the outermost field line which encloses 8.8 • 10 21 Mx 
flux. Note that both simulations have different domain 
sizes (49.152 x 49.152 x 16.384 Mm 3 in Figure [I] and 
73.728 x 73.728 x 9.216 Mm 3 in Figure^. We present 
a common subsection in both figures for Detter compar- 
ison. Figure [5] presents a 12 hour average from t = 12 to 
t = 24 hours, an intensity image at t = 24 hours for this 



simulation is presented in Figure |10| 

In both cases (Figure |4] sunspot without penumbra; 
Figure [5j sunspot with penumbra) the dominant fea- 
ture is a deep reaching outflow region surrounding the 
spot with radial mean velocities of about 0.5u rms (pan- 
els a). The corresponding large scale upflow (panels b) 
coincides with a large scale temperature enhancement of 
about 0.5T rms (panels d), indicating a thermal signature 
typical for a convective flow pattern. The amplitude of 
the radial flow as well as thermal signature is larger for 
the sunspot with penumbra. The same is also true for 
the fraction of the overturning mass flux found in the az- 
imuthal mean component (panels c) which is about 55% 
(80%) for the sunspot without (with) penumbra in 8 Mm 
depth. 

The most obvious difference between both sunspot 
models are present in photospheric layers along the indi- 
cated r = 0.1 level (see also Figure [6] for further detail): 
in absence of a penumbra we find a converging flow for 
R < 10 Mm and a diverging flow further out; in presence 
of a penumbra radial outflows are found everywhere in 
the photosphere. In absence of a penumbra the converg- 
ing flow shows a downward extension in the proximity 
of the sunspot, which disappears in the presences of a 
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Figure 6. Comparison of the radial flow amplitudes in the upper most 3 Mm of the domain. Panel a) corresponds to the simulation 
presented in Figure [4] panel b) to Figure [5] In contrast to Figure [4] and Figure [5] the velocity is not normalized by the rms velocity. Also 
note the different horizontal scales and resulting aspect ratio in both panels. The two dark lines indicate the t = 1 and r = 0.01 levels 
in both panels. In both cases the flow amplitudes of the outflows surrounding the sunspots peak above t = 1, typical flow velocities are 
about 600ms — 1 for the case without penumbra and 400ms -1 for the case with penumbra. The numbers indicate flow components that 
are discussed further in the text. 

the previous figures we did not normalize here the veloc- 
ity with the convective rms velocity. In the case of the 
sunspot without penumbra (panel a) we see outflows at 
photospheric levels for R > 11 Mm of about 600 ms -1 . 
These flows are a direct continuation of the subsurface 
flow cell identified above and they reach their peak am- 
plitude in the photosphere in between the r = 1 and 
t = 0.01 level. In the case of the sunspot with penumbra 
(panel b) the upward continuation of the deep reaching 
flow cell coincides with the penumbra and Evershed flow 
reaching an amplitude of about 4kms~ 1 . In addition 
we see at greater distance (R > 20 Mm) an additional 
superficial outflow pattern with amplitudes of about to 
400 ms -1 mostly above r = 1. This flow cell is neither 
related to the deep reaching flow component nor a con- 
tinuation of the Evershed flow. 

Figure [7] shows vertical profiles of this flow component 
for both sunspots together with the radial magnetic field 
and temperature. The largest flow amplitudes are found 
mostly above r = 1 and the velocity peaks in the region 
with the steepest increase of Br toward the overlying 
magnetic canopy. Since most of the stronger magnetic 
field is found in the overlying canopy, the convection pat- 
tern remains close to isotropic. This differs from the sit- 
uation found within the penumbra, where the magnetic 
field is strong en ough to cause substantia l anisotropy at 
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Figure 7. Vertical profiles of the radial magnetic field strength 
(black), radial flow velocity (blue) and temperature (red). Solid 
lines show averages in between R=13 and R=20 Mm in Figurc[6ji, 
dashed lines show averages in between R=23 and R=28 Mm m 
Figure [6b. The fastest outflow velocities are found about 200 km 
above the t = 1 level, where the radial magnetic field strength 
increases and forms a magnetic canopy overlying the plage region 
surrounding the sunspots. 

penumbra: radial outflows are found at all depth levels. 
As a consequence we do not find a situation in which a 
shallow Evershed flow is stacked on top of an inflow cell 
deeper down as it has been sugg ested by some helioseis- 
mic inversions (|Zhao et al.|2010[ ). 

Note that in i 'igure H] downfiows appears to be located 
within the sunspot, which is caused by the azimuthal av- 
eraging: the flow patterns presented in Figs. [2] and [3] 
clearly show that downflows are located at the periphery 
of the highly fragmented subsurface magnetic field struc- 
ture. The conspicuous temperature enhancement present 
in Figure [5]l) at the bottom boundary near R = 5 Mm 
is a consequence of our boundary condition, which keeps 
the magnetic field fixed within the inner most 8 Mm. 

Radial flows in the upper most 3 Mm of the domain 
for both cases are presented in Figure [6j In contrast to 



and below r = 1 |Kitiashvili et al.| (|2009|) ; IRempel et al 



( 2009a[ ); |Rempel| ( |201ip and leads to a radial flow veloc 
ity that peaks near r = 1. In that sense the flow shown 
in Figure [7]is better characterized as overshooting granu- 
lation that is outward deflected by the overlying inclined 
magnetic canopy. 

In addition strong inflows are present in layers higher 
than r = 0.001 regardless of the presence or absence of 
an penumbra, which are possibly related to the inverse 
Evershed flow ( Dialetis et al. 1985 ) . This flow component 



is very robust and present in all numerical sunspot simu- 
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Figure 8. Temporal evolution of large scale flows around the sunspot with penumbra. The quantities shown are similar to Figure [6] In 
contrast to Figure[(ib) where we presented a time average from 12-24 hours, we show here in panel a) a time average from 6-12 hours and in 
panel b) the average from 18-24 hours. Overall we see a trend of decreasing flow velocities in the region R > 20 Mm. While the photospheric 
flow speeds in panel a) are comparable to those found for the sunspot without penumbra (Figure pa), panel b) shows a reduction by about 
a factor of 2. In addition the separation of the flow component 4 (as defined in Figure |6| becomes more evident. 



lations we performed to date, however, due to the rather 
crude treatment of physics in those layers a more de- 
tailed study of this component would be required before 
stronger conclusions can be drawn (see also comments in 
Rempel et al.||2009a| |Rimpei1[20lT| ). 

Overall we identified in our simulations a total of 5 
flow components that contribute to horizontal flows in 
the proximity of sunspots. We highlighted the following 
flow components in Figure [6j (1) The fast photospheric 
Evershed outflow in the presence of a penumbra, (2) a 
converging fast photospheric inflow (with a weaker down- 
ward extension) in the absence of a penumbra, (3) a deep 
reaching outflow cell with average velocities of about 50% 
of the convective rms velocity, (4) a superficial outflow 
underneath the magnetic canopy, and (5) an inflow above 
t = 0.001. The flow components (1) and (2) are obvi- 
ously strongly affected by the presence or absence of an 
penumbra, while the flow components (3) - (5) are of 
independent origin. 

While the temporal evolution of the large scale flows 
surrounding the sunspot without penumbra (see Figure 
[2j right panels) shows an overall trend toward increas- 
ing outflow velocity, we find the opposite for the sunspot 
with penumbra. Figure [8] presents time averages of the 
radial flow velocity from 6-12 and 18-24 hours. The am- 
plitude of flows in the photosphere is declining by about 
a factor of two. While in earlier stages the flow compo- 
nent (4) showed some connection to the deeper reaching 
component (3), it is mostly disconnected toward the end 
of our simulation run. We return to this aspect in Sect. 
1531 



3.3. Origin of large scale subsurface flows 

The large degree of similarity between Figs. [4] and [5] 
strongly suggests an origin of the deep reaching flow com- 
ponent (3) that is independent from the Evershed flow. 
It is however non-trivial to clearly disentangle in the nu- 
merical simulations presented here the different processes 



that can give rise to large scale flows. There are primar- 
ily two effects: 1.) The sunspot imposes a geometric 
constraint on the surrounding convection pattern, which 
forces adjacent convection cells to align in a ring like pat- 
tern resulting in a mean flow. 2.) Reductions in the sur- 
face brightness lead to less downflows in the proximity of 
the spot. Since downflows provide low entropy material 
the result is an increase of the mean temperature driv- 
ing a large scale flow pattern through buoyancy. We see 
evidence for both processes at work. Figs. [2] and [3] show 
clearly the ring like arrangement of convection cells and 
also the strong reduction of downflow lanes within this 
region. To get an independent estimate of the strength 
of both processes we report here on a series of idealized 
simulations in a 49.152 Mm wide and 8.192 Mm deep 
domain. 

In the first setup we cut out a central cylinder with 
a radius of 8 Mm to independently investigate the geo- 
metric effect. To this end we artificially set all velocities 
to zero within this region and replace the stratification 
with the mean stratification found outside this cylinder. 
In the second setup we set the radiative cooling function 
to zero within a central cylinder of 8 Mm to investigate 
flows resulting from the blockage of heat flux. A combi- 
nation of geometrical and thermal effects can be realized 
by considering a cone instead of a cylinder. 

The resulting azimuthally averaged flows are presented 
in Figure [9j In the case of the pure geometric constraint 
(panel a) we see the development of a large scale flow 
cell, with flows converging toward the central cylinder 
for R < 12 Mm and diverging flows further out. Over- 
all the extent of the region with mean flows is about 
20 Mm, which is about 2.5 times the radius of the cen- 
tral cylinder. However, the amplitude of the outflow is 
with 10 — 15% of the convective rms velocity rather weak, 
larger values are found in the inflow cell in close prox- 
imity of the central cylinder. Heat flux blockage (panel 
b) leads to outflows at all depth levels in the proximity 
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Figure 9. Comparison of large scale flow patterns resulting from a combination of geometric constraints and blockage of heat flux in the 
photosphere. Panel a) shows the flow pattern developing around a cylinder with 8 Mm radius due to the imposed geometric constraint. 
Panel b) shows the flow patterns developing in response to switching off the radiative cooling in the region R < 8 Mm. Panels c) and d) 
show flows around cone shaped obstacles with different opening angles, which result from a combination of geometric and thermal effects. 



of central cylinder, even though further out we see the 
development of an inward directed flow cell. The mean 
flow speed is comparable to convective velocity in the up- 
per most 2 Mm but drops to values of less than 50% of 
the convective rms velocity in deeper layers. In this ex- 
periment we assumed complete blockage of the heat flux, 
which clearly overestimates the influence of a sunspot. A 
combination of both effects can be realized by placing an 
upward opening cone as obstacle into the center of the do- 
main (bottom panels) . The mean flow velocity increases 
with opening angle and values of about 50% of the con- 
vective rms velocity can be achieved with opening angle 
of 30 deg (panel d) . From these idealized experiments we 
can estimate that the bulk of the outflows results from 
blockage of heat flux, while geometric effects can have a 
contribution of up to 30%. 

The results from the idealized experiments show a good 
qualitative agreement in terms of amplitude and outflow 
pattern with the flows presented in Figs. H and [5] (pri- 
marily flow component 3). In the case oi the sunspot 
with penumbra (Figure [5]) the region with enhanced sub- 
surface temperature is as expected centered underneath 
the penumbra with reduced heat loss. Interestingly the 
region with enhanced temperature in (Figure El) is found 
in the same location, despite the absence ofa penum- 



bra. Furthermore the subsurface flow cells have in both 
cases about the same radial extent despite the fact that 
the sunspot with penumbra has almost twice the diam- 
eter. This is likely a consequence of the converging flow 
cell surrounding the sunspot without penumbra, which 
is shifting the location of diverging flow pattern further 
outward relative to the spot boundary. The converging 
flow cell is not present in the sunspot with penumbra. 

3.4. Are there bright rings around sunspots? 

Whether sunspots are surrounded by bright rings that 
account for some of the energy flux blocked by the lower 
brightness in umbra and penumbra has been extensively 
discussed in the literature. We refer here to [Rast et ahl 



(2001 1 and references therein for a detailed discussion of 
both theoretical and observational aspects of this prob- 
lem. On the theoretical side it has been argued that the 
thermal conductivity and heat capacity of the convec- 
tion zone are sufficiently large to absorb temporary ther- 
mal perturbations caused by sunspots without be coming 
visible as brightness change in the photosphere dSpruit 
1982a bl |Fbukal et al.|1983"l|Chiang fe Foukal|1985[|!jpruit, 
2000 ) . On the observational side the main challenge 
lies m the proper separation of this effect fro m Facu- 
lar br i ghten ing, which is a pure surface effect. Fowlei 
et al.| ( |1983[ ) found bright rings in the 0.1 — 0.3% range 
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Sure 10. Intensity image and magnetogram of the sunspot with penumbra for whic h w e present the subsurface flow structure in Figures 
and [8] long-term averages of the photospheric brightness are presented in Figure [TT] Displayed are intensity and magnetic field after 
running trie simulation for 24 hours. An animation is available in the online material. 



the more recent study by Rast et al. (20011 brightness 
enhancements of 0.5 — 1%, which are attributed to con- 
vective origin. 

We analyze here the numerical simulation in the 
73.728 x 74.728 x 9.216 Mm 3 domain, which has a more 
realistic penumbra and a moat region of at least 20 Mm 
surrounding the spot. A snapshot of t he inte nsity for a 
vertical ray is presented in Figure [TOj Figure [TT] shows 
azimuthally averaged intensity profiles for this sunspot 
model. We show three consecutive 6 hours temporal av- 
erages of the bolometric intensity for a vertical ray. Panel 
a) shows intensity profiles from the umbra into the moat 
region, panel b) fluctuations within the moat region on 
a different scale. We do not see evidence for a system- 
atic variation of intensity with radius in the moat region. 
For the individual 6 hour averages systematic trends are 
smaller than 0.1% and these trends are different for con- 
secutive 6 hour averages. The dark solid line shows the 
average of the intensity from 6 — 24 hours. Here we see 
a weak (0.05%) enhancement of intensity in the moat re- 
gion with a peak at about R = 30 Mm, however, given 
the substantial variation between the 6 hours averages it 
would require a longer simulation run to determine the 
robustness of this trend. The presented intensity profiles 
are based on azimuthal averages that take into account 
all magnetized areas in the moat region. Since our nu- 
merical resolution of 48 km in the horizontal direction 
in combination with grey radiative transfer is not suf- 
ficient for resolving Faculae, we do not have to correct 
for Facular brightening. Nevertheless, we repeated the 



analysis after cutting out magnetized regions. If we con- 
sider only regions with \B Z \ < 500 G the results remain 
essentially the same. If we lower the threshold to 250 G 
we see a bright ring with about 0.3% amplitude. Since 
most vertical magnetic field is concentrated in the darker 
downflow lanes and more unsigned vertical flux is present 
in the proximity of the sunspot the latter is an artifact 
of removing more dark downflow lanes from the average 
in the inner moat region compared to the outer. 

We also note that the total energy flux in our domain is 
reduced compared to quiet sun convection, since the con- 
vective energy flux crossing the bottom boundary is free 
to adjust (we specify the entropy in inflow regions, but 
the mass flux crossing the boundary is not constrained). 
In addition most of the energy flux that is blocked by 
the umbra of the sunspot does not enter the computa- 
tional domain in the first place due to our choice of the 
bottom boundary condition. Nevertheless, some of the 
energy flux entering the domain underneath the penum- 
bra has to be diverted away from the spot, which is, as 
we discussed above, the primary reason for the presence 
of the large scale subsurface flow cell surrounding the 
sunspot. The latter process is associated in the deeper 
layers with a temperature excess reaching about 50% of 
the convective rms temperature fluctuation (see Figure 
[51). The absence of bright rings at photospheric levels 
clearly shows that this temperature perturbation remains 
well hidden beneath the photosphere. 

4. INFLUENCE FROM BOTTOM BOUNDARY CONDITION 
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Figure 11. Temporal averages of the the azimuthally averaged intensity for the sunspot presented in Figure [To] Panel a) shows the 
Intensity profile from the center of the spot to the edge of the moat region, panel b) presents fluctuation around the mean intensity in the 
moat region. Dotted lines indicate the fluctuations present in 6 hour azimuthal averages. Solid lines are in addition averaged with a 2.5 Mm 
window in the radial direction to highlight potential trends of intensity with radius. The black line shows the average from 6 — 24 h. 

The simulation we reporte d on so far in S e ct. [3] as well 
as previous simulatio ns by Rempel et al. (2009b) and 
Rempel et al. ( 2009a ) were based on a bottom boundary 
that sets all velocities to zero in regions of strong field, 
but allows for convective motions crossing the boundary 
in regions with weak field. In the latter, velocity com- 
ponents are symmetric across the boundary in downflow 
regions, while the velocity is vertical (horizontal com- 
ponents antisymmetric) in upflow regions. In shallow 
domains with about 6 Mm depth this boundary condi- 
tion is a necessity to prevent a rather quick decay of a 
sunspot within a few hours. We demonstrate this here by 
a numerical experiment performed in a 6.144 Mm deep 
domain using an open bottom boundary, which imposes a 
symmetric boundary for all three mass flux components. 
In regions with \B Z \ < 2.5 kG we impose a constant 
total pressure at the bottom boundary, in regions with 
\B Z \ > 2.5 kG we extrapolate the pressure stratification. 
The latter is necessary in a shallow domain to prevent 
long lasting inflows that would destroy the sunspot on 
even sho rter time scales than reported here. 

Figure [12] shows the evolution of the sunspot in the 
6.144 Mrndeep domain. We restarted this simulation 
from a run which was previously evolved with a closed 
boundary condition in regions with \B Z \ > 2500 G, lead- 
ing to a near ly axisymmetric stable sunspot (top panels 
of Figure 12 ). After only a few hours we observe a decay 



of the spot that is primarily driven through interchange 
instabilities. After 3 hours we see the clear evidence of 
intrusions of field free plasma near the bottom boundary, 
which pushes buoyantly upward within the magnetized 
region and becomes visible in the photosphere through 
granulation patches showing up at the periphery of the 
umbra. After 6 hours these regions have expanded to 
patches of 5 — 10 Mm diameter and most of the umbra 
is dispersed along the periphery of the patches. 



In Figure [13] we contrast this result with a repetition 
of the same experiment in a 16.384 Mm deep domain. 
To this end we start from the snapshot at t = 8.8 hours 
shown in Figure [I] and evolve the simulation for 37.5 
hours with the open boundary condition. In contrast 
to the experiment performed in the 6.144 Mm deep do- 
main we impose a constant total pressure everywhere in 
the domain. Figure [13] presents 3 snapshots about 12.5 
hours apart from each other, the quantities shown are 
the same as in Figure [l] Suddenly opening the bot- 
tom boundary beneath the sunspot leads to an upflow 
strongly weakening the magnetic field, however, due to 
the deeper domain this upflow does not continue into the 
photosphere. After the initial weakening in the deep do- 
main the trunk of the sunspot fragments into smaller flux 
bundles. The photospheric signature of this process is an 
increased rate of deformation and fragmentation of the 
sunspot compared to the reference in Figure [T] with the 
closed boundary, but even at t = 46.3 hours (37.5 hours 
of evolution with open boundary condition) most of the 
flux remains concentrated together. The overall slower 
decay rate in the deeper domain compared to the 6 Mm 
deep domain discussed above is consistent with the in- 
crease of the convective time scale by a factor of about 
6 — 8 compared to the shallow domain. In Figure [14] 
we compare the flux remaining in the sunspots toward 
the end of our simulation runs. To evaluate the flux 
of the sunspots we compute the flux within the shown 
1.25 kG contour line, which is based on a smoothed 
magnetogram (Gaussian with FWHM of 960 km). In 
the case of the closed (open) boundary condition about 
86% (71%) of the initial flux at t = 8.8 hours (shown 
in panel a) is still concentrated together. At the same 
time the average field strength is dropping from 3300 G 
to 3040 G (2760 G), so that the actually remaining spot 
area is 93% (85%). Overall the decay rate is about a 
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Figure 12. Temporal evolution of a sunspot in a 6.144 Mm deep domain. The bottom boundary is symmetric in all three mass flux 
components. The sunspot is decaying on a time scale of about 3 — 6 hours due to interchange instabilities developing near the bottom 
of the domain. The field free plasma intruding near the bottom boundary leads to buoyant upflows that develop islands of granulation 
withing the umbra of the spot. After about 6 hours most of the magnetic field is dispersed into structure of the size of a few pores. 



factor of 2 larger for the open boundary condition. Av- 
eraged over the 37.5 hours of the simulated decay the 
corresponding flux loss rates are 0.8 • 10 21 Mx day -1 
(1.5 • 10 21 Mx day -1 ). These val ues are about a factor o f 
10 larger than those observed by Martinez Pillet ( 2002[ ), 
who found values around 0.6 — 1.44- 10^ u Mx da y -1 . For a 



substantially larger sunspot |Kubo et al.| ( |2007[ ) found val- 
ues of 0.6 — 0.8 • 10 21 Mx day -1 . A detailed comparison 
to observations is currently only of limited practicality 
since our simulated sunspots shown in Figure [l] and 13 
are essentially penumbra free "naked" sunspots and most 



reported sunspot decay rates refer to spots with fully de- 
veloped penumbrae. However, a direct comparison be- 
tween both simulations gives important clues about the 
influence of the bottom boundary condition. 

Since we start our simulations from a thermally relaxed 
HD simulation with magnetic field added to it, initially 
magnetic field and surrounding flow patterns are not nec- 
essarily dynamically consistent. In the case with closed 
boundary conditions (Figure [l]) the flow field is forced to 
adjust to the presence of the field, while in the case of 
open boundary conditions the magnetic field adjusts to 
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Figure 13. Temporal evolution of a sunspot in a 16.384 Mm deep domain using an open bottom boundary condition (symmetric in all 
three mass flux components). The simulation was restarted from the snapshot at t = 8.8 h shown in FigurelX] the elapsed time is relative to 
t = h to be comparable to Figure [l] The open boundary condition leads to a stronger deformation and fragmentation, compared to the 
reference r un s hown in Figure [T] but even at t = 46.3 hours the larger fr acti on of the flux remains concentrated within the main sunspot 
(see Figure [T4t . Compared to the 6 Mm deep domain presented in Figure (12) dynamical time scales near the bottom are about 6 — 8 times 
longer. An animation is available in the online material. 



the flow field, explaining the subs tant ial change in over- 
all shape of the sunspot in Figur e [T3| As a consequence 
the experiment shown in Figure |13| likely overestimates 
the decay rate due to this initial adjustment state. 

Overall our set of experiments demonstrates that the 
constraints imposed by the bottom boundary condition 
on the time evolution of sunspots can be relaxed by per- 
forming simulations in deeper (more computationally ex- 
pensive) domains. A 16 Mm deep domain is consistent 
with sunspot life times of a few days, which is still much 



shorter than the life time of most observed sunspots. Ex- 
trapolating this result a life time of a week or more should 
be reached in the 30 — 50 Mm depth range. The overall 
decay rate is likely also strongly dependent on the ini- 
tial state, which can be obtained only in a self-consistent 
manner through a flux emergence simulation describing 
the sunspot formation process. 

5. DISCUSSION 
5.1. Flow components around sunspots 




Figure 14. Influence of bottom boundary condition on sunspot decay. Panel a) shows the initial state at t = 8.8 h from which we started 
the simulation run with open boundary condition, panel b) shows the result for the closed boundary condition at t = 46.4 h, and panel c) 
for the open boundary condition at t = 46.3 h. Presented is a magnetogram at r = 1 together with a mask we used to compute the flux 
and area of the spot. The values are a) 9.1 ■ 10 21 Mx, b) 7.8 • 10 21 Mx, and c) 6.5 • 10 21 Mx for flux and a) 275 Mm 2 , b) 257 Mm 2 , and c) 
236 Mm 2 for the area. 



We find large scale outflows around sunspots as a ro- 
bust result from a scries of numerical simulations includ- 
ing idealized experiments. We identified 5 different flow 
components: 

1. The Evershed flow closely related to penumbral fine 
structure and driven by magneto convection in an 
inclined magnetic field. 

2. In the absence of a penumbra a converging flow 
caused by enhanced radiative cooling of granules 
adjacent to the umbra (flow extends with reduced 
amplitude to deeper layers). 

3. A deep reaching large scale flow that originates 
from a combination of geometric constraints and 
blockage of heat flux, with the latter providing the 
dominant contribution. The underlying convective 
flow morphology in the periphery of the sunspot is 
a ring-like arrangement of convection cells, which 
leads to large scale outflows reaching about 50% of 
the convective rms velocity over the depth range 
currently accessible through numerical simulations 
(down to 16 Mm depth). 

4. A shallow outflow component found in the plage re- 
gion surrounding the sunspots. This flow is mostly 
limited to the photosphere and is best characterized 
as overshooting granulation that is outward de- 
flected by the overlying inclined magnetic canopy. 

5. An inflow in levels higher than t = 0.001 located 
near the outer edge of the umbra or above the 
penumbra (if present). 

By comparing numerical simulations of sunspots with 
and without penumbra we identified an origin of the flow 
components (3) to (5) that is independent from the ex- 
istence of a penumbra and Evershed flow. Through a 
series of idealized experiments we demonstrated that the 
flow component (3) arises from a combination of geomet- 
ric constraints and blockage of heat flux, with the latter 
playing the dominant role. 



In models without penumbra we find a two cell flow 
pattern around the sunspot, consisting of a converging 
flow in the proximity of the spot (2) and a diverging flow 
further out (3), which merges in the photosphere with 
the flow pattern (4). The dominant flow component in 
this case is (3). 

In models with penumbra the flow pattern (2) is not 
present and replaced with (1). Flow pattern (1) merges 
continuously in the deeper layers with (3), resulting 
in outflows at all depth levels underneath the penum- 
bra. Flow pattern (4) is in this case somewhat detached 
and leads to a superficial flow component outside the 
sunspot. The amplitude of this flow component is de- 
clining throughout the simulation from about 600 ms _1 
to less than 300 ms _1 . 

The inflow component (5) in higher layers is present 
in both c ases and possibly rel ated to the inverse Ever- 
shed flow Dialetis et al. ( 1985 ) . This flow is very robust 
and present in all numerical sunspot simulations we per- 
formed to date, however, due to the rather crude treat- 
ment of physics in those layers a more detailed study of 
this component would be required before st ronger conclu- 
sions can be drawn ( see also comments in Rempel et al. 
M0^[Ftemiell[20lT| ). 

The Fvershed flow component (1) has been discussed 

with great detail in previous publications about MHD 

simulations (iHeinemann et al. 120071 iScharmer et al. 2008 
t> i „t „i'i'vwwiliiuu' „u 'i; J'i inA/wiiiu i' + „i 



Rem pel et ar.|2009b[|Kitiashvhi et al.|2009t|Kem~ 

2009a Rempel 2011) and is well studied in observ 



:mpcl 2011) 
see e.g. |Solarhrij|2UU3| |Thomas fc Weiss 



et a 
ser vat ions 
20041), we focus 



the following discussion primarily on the flow compo- 
nents (2) to (4). 

The inflow component (2) has be en observed in pores 
as well as spots without penumbra (Wang fc Zirin||1992 
et al.| [T999; Vargas Dommgu ez et al.||2010p anc 
i identified in 3D MHU simulations of pores 



Sobotka 



has been 



(Cameron et al. 2007), where it plays a crucial role in 
maintaining the magne tic structure as well as forming it 
< |Kitiashvili et al.||2010| ). 

The dominant subsurface outflow component (3) is es- 
sentially a convective flow that is modified due to a com- 
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bination of geometric constraints together with heat flux 
blockage caused by the presence of the sunspot. As- 
sociated with this flow is a ring of temperature excess 
around the sunspot, which is similar to the flow pattern 
an d related thermal perturbations that were suggested 
by Meyer et al. ( 1974 ) , although we see this flow on scales 
smaller than supergranulation. Despite a thermal sig- 
nature in deeper layers with amplitudes comparable to 
those of convective rms temperature fluctuations we do 
not see compelling evidence for a bright ring surrounding 
the sunspot in the photosphere. 

Apart from being a shallow photospheric moat flow 
component, the flow pattern (4) could be possibly also 
rela ted to observed flow s along the magnetic canopy (see 

and references therein) 



e.g 



Rezaci ct al. 2006 



Since 

we see this effect in our simulations regardless of the 
presence or absence of a penumbra this component is 
not an extension of the Evershed flow in our simulations. 

The flow pattern we find around the sunspot without 
penumbra (converging flow in close proximity, diverging 
flow further out) is similar to flo ws found in axisymmet- 
ric su ns pot simulations suc h as IHurlburt fc Ruc klidgc 
(pOO] ), [Botha et al| ( [2006] ), and |Botha et al.| (|2008). 
It has been speculated that in the case of a sunspot 
with penumbra only a shallow outflow (Evershed flow) is 
added on top of this flow pattern, whi l e a co nverging flow 
remains in deeper layers Zhao et ah] (2010). We do not 
find any evidence for such a layered now structure in our 
simulations. In the case of a sunspot with penumbra the 
inflow cell adjacent to the spot disappears and the spot 
is surrounded by outflows at all depth levels (a combina- 
tion of the flow patterns (1) and (3)). This is consistent 
with the underlying driving mechanisms of these flows: 
In the absence of a penumbra granules in the proximity of 
the umbra have enhanced radiative losses, creating cool 
downflows that extend to the bottom of the domain and 
maintain in the proximity of the sunspot a converging 
convective flow pattern. In the presence of a penumbra 
the region with enhanced radiative loss is replaced by 
a region with reduced radiative loss, which leads to an 
overall reduction of cool downflows. In addition the fast 
near surface Evershed flow leads to a preferred draining 
of these downflows toward the outer edge of the penum- 
bra. The consequence is a temperature enhancement in 
the subsurface layers driving a broad upflow underneath 
the penumbra that results in outflows at all depth levels. 
Th is flow patterns is consiste nt with a recent inversion 
by Gizon et al. (2009 2010b I that reports on outflows 
in the upper mos t 4.5 Mm. A recent result by |Feather-] 
stone et al. ( 2011 ) indicates two contributions to the large 
scale flows surrounding sunspots: a superficial near pho- 
tospheric and a deeper reaching component that peaks 
at about 5 Mm depth and essentially disappears in depth 
greater than 9 Mm. 

Our sunspot models are set up in way that they yield 
almost quasi-stationary solutions, which implies that 
they are likely best compared to very stable, almost cir- 
cular sunspots with little evolution over time scales of a 
few days. It is possible that the subsurface flow structure 
is much more complicated around more complex active 
regions with rapid evolution. 

5.2. Connection to observed moat flows 



In observations moat flows are typically defined as out- 
flows that are present from the outer edge of a sunspot 
(penumbra) to about 2 sun spot radii. Typical flow yeloc- 
ities are around 500 ms _1 ( Brickhouse &: Labonte | ^988 - 
Sobotka fc Roudier|[2007] |Balthasar &: Muglach||2010[ ). 



Over the past years it has been discussed in the litera- 
ture whether Evershed and moat flow are connected or 
of ind epen dent origin. fSainz Dalda fc Martinez Pillet| 
( 2005 1 and Cabrera Solana et al.: (2006 ) found examples 
of moving magnetic features in the mo at region that show 
a connection to penumbral filaments. |Vargas Dommguez 



et al. ( 2008 ) studied several examples of complex active 



regions and found a very close connection between the 
presence of penumbrae and moat flows. Furthermore, 
moat flows were only found near penumbrae perpendic- 
ular to the sunspot border and absent near tangential 
penumbrae. Moat flows were also absent when there was 
a polarity inversion lin e in the proximity of the sunspot. 



Zuccarello et al. ( 2009 ) reported on moving magnetic fea- 
tures diverging from a sunspot without penumbra, indi- 
cating an origin of horizontal outflows independent from 
a penumbra and Evershed flow. 

In the simulations presented here we find sunspots sur- 
rounded by outflows in the photosphere with velocities in 
the 300 — 600 ms _1 range that extend to about 2 sunspot 
radii, regardless of the presence or absence of a penumbra 
and the associated Evershed flow. 

In the absence of a penumbra we find in close proximity 
of the spot a ring of inflowin g plasma, which has been also 



observed ar o und pores by | Wang fc Zirin| (fl9 92|) ; |Sobotka| 



M ^l j|992 

et al. ( 1999|) ; |Vargas Dommguez etall (pOlo l ) . In addition 
Sobotka et al.| ( |l999[ ); | Vargas Dommguez et al. (20101 
observed a diverging flow further away that is related 
to a ring-like arrangement of "centers of positive diver- 
gence" around the pore, which is similar to the preferred 
ring-like arrangement of convection cells we find around 
our simulated sunspots over a depth range of more than 
10 Mm. In that sense observations likely see the "tip 
of the iceberg" of this flow structure. Large scale out- 
flows are also observed around " naked" (penumbra free) 
sunspots ( Zuccarello et al.|2009 ), which we would expect 
as this process is largely scale invariant. 

Interestingly, our sunspot with penumbra has an on 
average a weaker outflow at photospheric levels com- 
pared to the penumbra free sunspot. In addition the 
outflow results primarily from the very superficial flow 
component (4) and we observe a trend of declining flow 
speed throughout the simulation, which is opposite to 
trend we find in the simulation with the penumbra free 
sunspot. At later stages (> 12 hours) of the simula- 
tion with penumbra we find a rather sharp decline of the 
Evershed flow and subsurface flow pattern (3) at about 
R = 22 Mm, which results in a downflow that collects 
most of the horizontal mass flux present. This down- 
flow is in part driven by cool material deposited there by 
the Evershed flow (cool plasma forming in the penum- 
bra is transported preferentially outward), and in part a 
consequence of opposite polarity magnetic flux accumu- 
lating near the outer edge of the penumbra that diverts 
horizontal flows downward. The lower temperature is ev- 
ident from Figure[5]jd), the opposite polarity flux causing 
a p olarity inversion line around the sunspot from Figure 
|10| The absence of moat flows in a region with a polar- 
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ity inversion ne ar the outer edge of the penum bra was 
also observed by Vargas Dommguez et al. (2008). If this 
shielding effect would be weaker it is conceivable that the 
Evershed flow component (1) and deep flow component 
(3) have also larger contributions further out, leading 
to a deeper reaching and faster moat flow around the 
sunspot with penumbra as it is partially indicated dur- 
ing earlier stages of this simulation (see Figure[8^). How- 
ever, it would require additional numerical experiments 
to clearly quantify the role of a magnetic inversion line 
in our simulations. 

While we find large scale outflows around sunspots as 
a robust result regardless of the presence or absence of 
a penumbra, it is also possible that we are still missing 
some processes contributing to larger scale flows, in par- 
ticular on scales of supergranulation that are still not 
well captured by the extent of our simulation domain. 
Independent from that the simulations presented here 
clearly show that there are likely different contributions 
to large scale flows like the moat flow, and some of these 
contributions are independent from the Evershed flow 
(components 3 and 4). Based on our findings we conjec- 
ture that there is likely no clear "yes" or "no" answer 
to the question of whether moat and Evershed flow are 
related. To some extent this connection is also subject to 
interpretatio n and the adopted definition of " moat flow" . 
For example Vargas Dommguez et al. (2010) did not re- 
fer to the outflow they observed around pores as moat 
flow, while we see strong indications that those flows have 
likely a much deeper reaching structure and are in part 
related to blockage and modification of convective energy 
transport around sunspots. 

Due to the substantial differences in the depth extent of 
the individual outflow components we identified (compo- 
nents (1) and (4) are superficial, while (3) is deep reach- 
ing), a clarification of the Evershed- moat flow relation- 
ship has to consider the deeper reaching subsurface flow 
structure throu gh helioseismology (see, fo r example, re- 
cent result s by Gizon et al. 2009 2010b[ Feathcrstonc 
et al.||2011| ). 



5.3. Robustness of flow structure 

The large scale flows we discussed in previous sections 
reach scales comparable to the computational domain 
size. It is therefore essential to discuss to which degree 
the flow structure could be influenced by the overall sim- 
ulation setup including boundary conditions. A strong 
dependence on boundary conditions has been found for 
examp l e in t he 2D axisymmetric simulations of Botha 



et al. (20081. We have shown already that the mag- 
netic boundary conditions strongly influence the proper- 
ties of the magnetic sunspot structure. While the bot- 
tom boundary primarily affects the long-term stability 
and evolution of the sunspot, the top boundary condi- 
tion sets the overall extent of the sunspot penumbra in- 
cluding the associated Evershed flow. With regard to 
large scale flows the key question is whether these flows 
are a response to the magnetic structure present in the 
simulation and are driven by resolved physical processes 
within the computational domain or whether they are in 
addition substantially influenced by the hydrodynamical 
boundary conditions. 

Our top boundary condition is located about 700 km 
above the quiet sun r = 1 level and therefore a density 



contrast of about three orders of magnitude away from 
the photosphere. Although flow velocities can be sub- 
stantial (larger than 10 km s _1 ) at this boundary, the 
associated mass and momentum fluxes are negligible and 
no substantial feedback on flows in the photosphere and 
below takes place. The only flow component potentially 
affected by this boundary condition is the component 
(5). Both, the Evershed flow (1) and inflow in the ab- 
sence of penumbra (2) arc driven by resolved physical 
processes within the computational domain. This is also 
the case for the flow component (4), which is in addition 
less dependent on the magnetic boundary condition than 
components (1) and (2). 

The largest scale flow we find in our simulations is the 
deep reaching flow component (3). In all of the simula- 
tions discussed this flow extends to the bottom bound- 
ary and at least in the domains that are only 50 Mm 
wide the radial extent of this flow becomes comparable 
to the horizontal domain size. The potential influence of 
the horizontal and bottom boundary conditions on this 
flow could be investigated by either changing the bound- 
ary conditions or by comparing simulations of similar 
sunspots in differently sized domains. We have done here 
the latter. Figures[4land[5]compare flow systems in 2 sim- 
ulations with quite (different domain size, only a common 
subsection is shown. In Figure [4 the simulation domain 
extends to about 15.5 Mm depth, i.e. the mass in the 
part not shown exceeds the mass in the part shown by 
about a factor of 10! Nevertheless the flow structure of 
component (3) is very similar in both cases, and in par- 
ticular the absence of a converging return flow near the 
bottom boundary is a robust result (such a flow was only 
present as a transient during early evolution stages shown 
in Figure [2| . Similarly the horizontal extent of the do- 
main does not influence the flow structure substantially. 
Even though the domain extends in Figure [5] to a radial 
position of 37 Mm, the radial extent of the large scale 
flow is similar to that found in Figure [4] We have done 
several additional simulations we did not report here, in 
particular a simulation of the sunspot shown in Figure [4] 
in an only 8 Mm deep domain and a series of simulations 
with different grid resolution of the sunspot shown in Fig- 
ure [5] in a domain only 49.152 Mm wide and 6.144 Mm 
deep. In all cases we find with regard to the axisymmet- 
ric flow components similar results when we compare the 
overlapping parts of the simulation domains. Additional 
evidence that the horizontal boundary conditions do not 
matter too much for this flow co mponent comes f rom th e 



"double sunsp ot" simulation of Rempel et al. (2009a); 



Rempel (2011). Due to the horizontal periodicity this 



simulation assumes sunspots of alternating polarity in 
the x and same polarity in the y direction. While this 
strongly influences the structure of the penumbra and 
Evershed flow, the subsurface flow pattern corresp onding 



to flow co mponent (3) remains mostly unaffected (Rem 



pei][20lT] see Figure 20) 



A related concern is the temporal evolution of large 
scale flows as our simulations address at this point mostly 
rather short lived sunspots. This is not a major concern 
for the flow components (1), (2) and (5) which are es- 
tablished within a few hours of simulation time and did 
not show any substantial variation over the time frame 
of 1-2 days we covered. The flow components (3) and (4) 
show a temporal evolution which is presented in Figures 
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tand [8] The rather strong variation of component (4) 
own in Figure 8] is limited to this particular sunspot 
and likely relatedto the magnetic fiel d ev olution in the 
moat region as discussed in Section |5.2| We did not 
observe a substantial variation of this how component 
around the sunspot without penumbra. Since the flow 
component (3) essentially extends from the photosphere 
to the bottom boundary of our domain it covers regions 
with a substantial variation of the intrinsic convective 
time scales. The quantity H p /v rms varies from about a 
minute in the photosphere to 6 hours in 15.5 Mm depth. 
The effective life time of convection patterns at the re- 
spective depths is about a factor of 5-10 larger, leading to 
time scales of 5-10 minutes in the photosphere (life time 
of granules), 3-6 hours in 5.5 Mm depth (approximate 
life time of the sunspot in Figure 12 1, and 30-60 hours in 
15.5 Mm d epth (approximate life time of the sunspot in 
Figure 13 ). From this we can conclude that the the large 



scale flow component (3) can be considered mostly con- 
verged in the upper most 8 Mm of our simulation domain 
(at least 5-10 convection pattern life times). In addition 
we did not see any evidence that the temporal evolution 
changes the large scale flow topology, it mostly affects 
the average flow amplitudes (see Figure [2]). As far as it 
concerns the flow component (3) the current simulations 
likely underestimate the flow velocity in the deeper parts 
of the domain. 

Related to the overall temporal evolution is also the 
initial state of our simulations and its influence on the 
subsequent sunspot evolution as discussed in Section |4j 
In our current setup large scale flows are the response to a 
magnetic obstacle in the convection zone that is initially 
not consistent with the flow fields present. However, the 
system evolves into a state in which both become con- 
sistent by evolving the magnetic field and flow structure. 
Ultimately this shortcoming will be resolved through ac- 
tive re gion scale flux emerge nce simulations. Fist results 
such as Cheung et al. ( 2010 ) and work in progress lead to 
results that are in terms of magnetic field and flow struc- 
ture comparable to the penumbra free sunspot discussed 
here. 

5.4. Subsurface field structure of sunspots 

We investigated the influence from the bottom bound- 
ary on the overall stability and life time of our simu- 
lated sunspots by comparing simulations in 6.144 and 
16.384 Mm deep domains. In the shallow domain a closed 
boundary condition in strong field regions is required to 
prevent an unrealistically fast decay of the sunspot. Us- 
ing the same boundary condition in the 16 Mm deep 
domain adds several degrees of freedom to the long term 
evolution of the sunspot allowing for flux separation and 
decay driven through convective motions in more than 
10 Mm depth beneath the solar surface. In the photo- 
sphere these flux separation events are accompanied with 
light bridge formation. Changing the bottom boundary 
condition to an open boundary in strong field regions 
enhances the deformation and decay rate of the sunspot, 
but unlike shallow domains, the sunspot stays coherent 
for at least 24 hours in time. It is very likely that the 
time scales of sunspot decay will become consistent with 
observed sunspot life times in deeper computational do- 
mains regardless of the bottom boundary condition used. 
We note that these conclusions are currently based on 



numerical simulations of sunspots without an extended 
penumbra. It has been suggested that the presence of 
an extended penumbra can lead in the upper most few 
Mm to an additional sta bilization of a sun spot against 
interchange instabilities (Meyer et al. 19771. Note that 
our simulation with penumbra was performed only in a 
9 Mm deep domain (due to the substantially higher grid 
resolution) and used again a closed boundary in strong 
field regions. 

The "anchoring problem" of sun spots (see for exam- 
ple the correspondi ng chapters in |Gizon et al.| [2010a; 



Moradi et al. 



2010 



and further references theremj is 
closely related to the subsurface structure of sunspots. 
A sufficiently deep anchoring is required to explain the 
relatively long life times of sunspots compared to the 
convective turnover as well as Alfvenic travel time scales 
found in the upper convection z one. The two solutions 



discussed are the cluster model (Parker 19791, in which 
a convergent subsurface flow opposes sunspot decay, and 
a sufficiently deep anchoring close to the base of the con- 
vection zone where convective turnover time scales (de- 
fined through H p /v rms ) are of the order of 1-2 weeks or 
even longer (overshoot region). Note that this is about 
a factor of 30 longer than the longest time scales in our 
16 Mm deep simulation domain, in addition we found 
that the effective life time of the magnetic structure is 
about 5 — 10 H p /v rms . 

With regard to the cluster model we do not see clear ev- 
idence that converging flows in the proximity of sunspots 
are a viable solution - at least not for a sunspot with an 
extended penumbra where we find outflows at all depth 
levels. Although a weaker converging flow is present in 
our sunspot model without penumbra, we do not see a 
clear indication that it stabilizes the sun spot sufficiently 
against decay (see the case in Figure 13 1. 

With regard to "deep anchoring" itnas been pointed 
out that this might be inconsistent with the photospheric 
motion of sunspots after the flux emergence process. Ris- 
ing magnetic flux tubes show typically rather long wave- 
lengths such as m = 1 or m = 2 modes (c.f. |Fan et al 



1993} [19941 IMoreno-Insertis et al.|[l994| |Schussler et al. 



1994| |Caligari et al. 1995p . Due to magnetic tension 



forces the opposite polarities of the active region continue 
to separate after the emergence at a rate that is inconsis- 
tent with observations. Because of this it has been sug- 
gested that sunspots become dynamically disconnected 
from their magnetic roots (at the base of the convection 



(Fan et al.||l994 


Moreno-Insertis et al. 


1995 


Schrijver 


& Title||l999| |Sc 


tnissler & Rempel||2005p. Here the term 



in which the magnetic field strength drops over a certain 
height range below equipartition field strength so that 
magnetic field becomes passive with respect to convec- 
tive motions. While the dynamical disconnection process 
alleviates the drift problem mentioned above, there are 
currently 2 unresolved issues: 1. Is it possibl e to dynam- 
ically di s conne ct a flux tube as suggested by Schiisslcr & 



Rempel| ( |2005[ ) without destroying the coherence of the 
sunspot in the photosphere? 2. Is the resulting shal- 
low (possibly less than 10 Mm deep) sunspot sufficiently 
stable? 

The models presented here give some clues to answer 
these questions. In the simulation with open bottom 



Subsurface magnetic field and flow 



19 



boundary condition in a 16 Mm domain we observe ini- 
tially (after switching from a closed to an open boundary) 
an inflow into the sunspot with an amplitude comparable 
to the typical convective rms velocity at the depth of the 
boundary condition, i.e. a velocity of a few 100 ms -1 . 
The latter is very similar t o the assumptions made by 
Schussler & Rempel (2005). The inflow leads to a de- 
struction of the coherent trunk of the spot in the lower 
part of the domain while the field in the photosphere 
shows only a moderate enhancement of decay compared 
to the closed boundary reference case. While this is in 
pri ncipal support of the discon nection scenario discussed 
by Schussler & Rempel (2005), our results point also to- 
ward the necessity of a sufficiently deep disconnection to 
ensure coherence of the photospheric parts of the sunspot 
over time scales of several days to weeks. In our cur- 
rently 16 Mm deep domain we observe a slow decay of the 
sunspot on time scales of days. Extrapolating this result, 
life times of about a week should be reached in a depth of 
about 30 — 50 Mm. In addition it is conceivable that the 
subsurface convection pattern is substantially altered as 
consequence of the flux emergence process forming the 
sunspots in the first place. This could lead to systematic 
difference compared to the simulations presented here, in 
which we added the magnetic field to a pre-existing con- 
vection pattern. This problem will be ultimately resolved 
once realisti c flux emergence simulations such as [Cheung 
et al.| ( |2010[ ) will become available in sufficiently deep do- 
mains, perhaps coupled with mod els of flux e mergence in 
the deep convection zone such as |Fan| ( |2008[ ) . 

With regard to the subsurface structure we point out 
that the numerical sim ulations presented her e and espe- 
cially those described in |Rempel et al. (2009a) are clearly 
biased toward the monolithic picture due to the mono- 
lithic initial state as well as bottom boundary condition 
(at least if the closed boundary is used). Nevertheless, we 
find in the deeper domains considered here a substantial 
fragmentation of the subsurface field caused by the inter- 
action with the surrounding convection. Almost all frag- 
mentation events (intrusions of field free plasma) present 
several Mm beneath the surface lead on time scales of a 
few hours to a day to the formation of light bridges or flux 
separation in the photosphere. These results allow for 
the conclusion that the subsurface field of sunspots very 
likely shows significant fragmentation, but most of these 
subsurface fragmentations (at least those present in the 
upper most 10 Mm) do not remain hidden and become 
visible in the photosphere on a rather short time scale, 
i.e. the photospheric appearance of a sunspot should tell 
us a lot about its subsurface structure: Sunspots that 
are very stable and do not show light bridges are more 
monolithic than sunspots with light bridges and signs of 
flux separation. 

6. CONCLUSIONS 

We presented simulations of sunspots in up to 16.384 
Mm deep and up to 73.728 Mm wide domains covering 
time evolution of up to 2 days. Through variations of the 
magnetic top boundary condition we simulated sunspots 
with and without penumbrae while having comparable 
magnetic flux, allowing for a direct side by side compar- 
ison of their properties with regard to large scale flows. 
We identified 5 different flow compo nent s in our simula- 
tions, which we discussed in Section 5.1 and highlighted 



possible connec tion s to the Evershed - moat flow connec- 
tion in Section IB~2l Overall we identified in our simula- 
tions 2 outflow components (one shallow and one deep 
reaching) that give a contribution to large scale flows 
around sunspots independent from the Evershed flow. 
In the case of a sunspot with penumbra we find outflows 
at all depths covered by the numerical simulation. Re- 
solving the Evershed moat flow connection requires in 
addition to the detailed study of photospheric flows also 
the study of deeper reaching flow structures that are ac- 
cessible through helioseismology. 

The long-term stability and evolution of the simulated 
sunspots is mostly governed by the longest time scale of 
convective motions found near the bottom of the simula- 
tion domain. If the magnetic field is not artificially con- 
strained at the bottom boundary, sunspots strongly de- 
form and decay on a time scale of the order of 10 H p /v rms 
(evaluated near the bottom of the domain). This trans- 
lates to a few hours in 6 Mm depth, about 1-2 days in 
16 Mm depth, and 10 days in 50 Mm depth (assuming the 
trend indicated in the numerical simulations shown here 
is valid over a broader depth range). We did not find 
a clear indication that the large scale flows developing 
around our simulated sunspots play a key role in main- 
taining their coherence, at least not beyond the intrinsic 
convective time scales present. In that sense our results 
point more toward "deep anchoring" as an explanation 
for the observed sunspot life times. Further clarification 
of this matter will require simulations in deeper domain 
that also include the flux emergence and sunspot for- 
mation process in order to start from a more consistent 
initial state. 
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